clear all;
clc;

Re = 6378.145e+3;

t0 = 0;
tf = 24*3600;
tspan = [t0, tf];

r0 = 1e+6 * [-0.24322284327422; -3.79199144866264; 6.32444994784559];
v0 = 1e+3 * [-1.64209188480826; 6.17220916529869; 3.63756062689151];
x0 = [r0; v0];

option = odeset('RelTol', 1e-6);
[t_orbit, x_orbit] = ode45('dif_orbit_dynamics', tspan, x0, option);
%[t_orbit, x_orbit] = runge('dif_orbit_dynamics',tspan, x0, 30);
%%
figure;
subplot(3,1,1);
plot(t_orbit, x_orbit(:,1));
xlabel('fly-time(s)');
ylabel('X(m)');

subplot(3,1,2);
plot(t_orbit, x_orbit(:,2));
xlabel('fly-time(s)');
ylabel('Y(m)');

subplot(3,1,3);
plot(t_orbit, x_orbit(:,3));
xlabel('fly-time(s)');
ylabel('Z(m)');
%%
figure;
subplot(3,1,1);
plot(t_orbit, x_orbit(:,4));
xlabel('fly-time(s)');
ylabel('Vx(m)');

subplot(3,1,2);
plot(t_orbit, x_orbit(:,5));
xlabel('fly-time(s)');
ylabel('Vy(m)');

subplot(3,1,3);
plot(t_orbit, x_orbit(:,6));
xlabel('fly-time(s)');
ylabel('Vz(m)');
%%
figure;
plot3(x_orbit(:,1), x_orbit(:,2), x_orbit(:,3),'r');
xlabel('X(m)');
ylabel('Y(m)');
zlabel('Z(m)');

%%
grid on;
hold on;
[Xe, Ye, Ze] = sphere(24);
Xe = Re * Xe;
Ye = Re * Ye;
Ze = Re * Ze;
surf(Xe, Ye,Ze);
axis equal;

%%
longitude_beijing = (116+28/60)*pi/180;
latitude_beijing = (39+48/60)*pi/180;
X_beijing = cos(latitude_beijing)*cos(longitude_beijing)*Re;
Y_beijing = cos(latitude_beijing)*sin(longitude_beijing)*Re;
Z_beijing = sin(latitude_beijing)*Re;
plot3(X_beijing,Y_beijing,Z_beijing,'pr');
